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Weak interactions of solitary waves in the generalized nonlinear Schrodinger equations are stud- 
ied. It is first shown that these interactions exhibit similar fractal dependence on initial conditions 
for different nonlinearities. Then by using the Karpman-Solov'ev method, a universal system of 
dynamical equations is derived for the velocities, amplitudes, positions and phases of interacting 
solitary waves. These dynamical equations contain a single parameter, which accounts for the dif- 
ferent forms of nonlinearity. When this parameter is zero, these dynamical equations are integrable, 
and the exact analytical solutions are derived. When this parameter is non-zero, the dynamical 
equations exhibit fractal structures which match those in the original wave equations both qualita- 
tively and quantitatively. Thus the universal nature of fractal structures in the weak interaction of 
solitary waves is analytically established. The origin of these fractal structures is also explored. It is 
shown that these structures bifurcate from the initial conditions where the solutions of the integrable 
dynamical equations develop finite-time singularities. Based on this observation, an analytical crite- 
rion for the existence and locations of fractal structures is obtained. Lastly, these analytical results 
are applied to the generalized nonlinear Schrodinger equations with various nonlinearities such as 
the saturable nonlinearity, and predictions on their weak interactions of solitary waves are made. 
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I. INTRODUCTION 



Solitary wave interactions are a fascinating and im- 
portant phenomenon for both physical and mathemati- 
cal reasons. Physically, such interactions have arisen in 
a wide array of disciplines such as water waves JM, op- 
tics @, S 13, H, H, 01, and Josephson junctions [§]. For 
instance, in soliton-based fiber communication systems, 
optical pulses traveling in different frequency channels 
pass through each other, giving rise to collisions (strong 
interactions) of solitary waves. In the same frequency 
channel, neighboring optical pulses interfere with each 
other through overlapping tails, giving rise to weak in- 
teractions of solitary waves. Motivated by these physical 
applications, solitary wave interactions has been stud- 
ied extensively in both the mathematical and physical 
communities. If the system is integrable, collisions of 
solitons are elastic [l| , and weak interactions of solitons 
exhibit interesting yet simple behaviors [2|, 13, LLC], [ill, LL2| ■ 
However, in non-integrable systems, solitary wave inter- 
actions can be far more complex. The first sign of this 
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complexity was reported by Ablowitz, et al. [13| for kink 
and anti-kink collisions in the cf> 4 model where, inside the 
trapping interval, a reflection window was found. Later 
extensive numerical studies on this model by Campbell, 
et al. [3, HE OH, 03 revealed that in fact, sequences of 
two- and more-bounce reflection windows exist, and the 
physical mechanism for these refection windows is a res- 
onant energy transfer between the translational motion 
and internal modes of kinks/antikinks. Anninos, et al. 
[HI pointed out further that there is a fractal structure 
in kink-antikink collisions. Using a collective-coordinate 
(i.e., variational) approach, they derived a set of fourth- 
order ordinary differential equations (ODEs) for these 
collisions, and these ODEs exhibit qualitatively similar 
fractal structures as in the </> 4 model (a comprehensive 
review on kink-antikink collisions in 4 -type equations 
can be found in [19(). These complex dynamics turn out 
to be not restricted to kink-antikink collisions. Indeed, 
similar phenomena have been reported on kink-defect col- 
lisions in the sine-Gordon and <fi 4 models [2(| [HI 
as well as vector-soliton collisions in the coupled nonlin- 
ear Schrodinger (NLS) equations [23l. l24l [2a|. Further- 
more, fractal scattering has also been reported on weak 
interactions of breathers in a weakly discrete sine-Gordon 
equation [2rj| and weak interactions of solitary waves in 
a weakly discrete NLS equation [27| . Recently, Good- 
man and Haberman [28|, [H, L13] provided a deep analy- 
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sis on the collective-coordinate models (ODEs) for kink- 
antikink collisions in the <fi 4 model 18j], kink-defect col- 
lisions in the sine-Gordon model [20 1, and vector-soliton 
collisions in the coupled NLS equations [25|, [3lJ , using so- 
phisticated dynamical system techniques. They derived 
analytical formulas for the locations of reflection- window 
sequences, which agree qualitatively with numerical re- 
sults on the original partial differential equations (PDEs). 
Their results shed much light on the origins of these win- 
dow sequences and fractal structures, especially from a 
mathematical point of view. 

Despite the above progress on solitary wave interac- 
tions, our understanding on these phenomena is far from 
satisfactory. On the collision of solitary waves, the anal- 
ysis done so far were all based on approximate collective- 
coordinate approaches, hence the reduced ODE models 
can only provide qualitative results at best. Many fea- 
tures reported in the ODE models can not be seen in 
the PDE simulations, thus it is not possible to make a 
reliable prediction on the collision dynamics based on 
those ODE models and their analysis. In addition, the 
ODE models obtained from the collective-coordinate ap- 
proaches not only are complicated, but also differ sig- 
nificantly from one PDE system to another. This forced 
previous researchers to analyze each PDE and its reduced 
ODE systems on an individual basis, which prevents an 
overall understanding on collision processes of solitary 
waves. On the weak interaction of solitary waves, the sit- 
uation is even less satisfactory. The fractal nature of this 
weak interaction was reported only for systems which are 
weakly perturbed integrable systems (sine-Gordon and 
NLS equations, to be more specific) [2(| [27j . It is not 
known yet whether similar phenomena arise in strongly 
non-integrable equations. More seriously, the previous 
work on this subject is largely numerical. No analysis has 
been attempted yet (not even the approximate collective- 
coordinate studies) . Thus an analytical understanding on 
weak interactions of solitary waves is a completely open 
question. 

In this paper, we study weak interactions of solitary 
waves in a whole class of generalized NLS equations (with 
arbitrary nonlinearities) both analytically and numeri- 
cally. These generalized NLS equations are not weak 
perturbations of the NLS equation in general. First we 
show by direct PDE simulations that these weak inter- 
actions for different nonlinearities exhibit similar fractal 
structures on initial parameters of solitary waves. This 
establishes that fractal scattering is a common feature of 
weak interactions in this class of generalized NLS equa- 
tions. Next, we rigorously derive a universal system of 
dynamical equations (ODEs) for the velocities, ampli- 
tudes, positions and phases of interacting solitary waves 
in this class of PDEs by the Karpman-Solov'ev method. 
This universal ODE system is remarkably simple, and it 
contains only a single parameter which depends on the in- 
dividual PDEs (after variable rescalings) . When this pa- 
rameter is zero, these dynamical equations are integrable, 
and their exact analytical solutions are derived. When 



this parameter is non-zero, the dynamical equations are 
found to exhibit fractal structures for a wide range of ini- 
tial conditions. These fractal structures match those in 
the original PDEs both qualitatively and quantitatively, 
thus the universal nature of fractal scattering in the weak 
interaction of solitary waves is analytically established. 
We further explore the origin of these fractal structures. 
Our numerical studies on the ODE system show that 
these fractal structures bifurcate from the initial con- 
ditions where the solutions of the integrable dynamical 
equations develop finite-time singularities. Based on this 
observation, we present an analytical criterion for the 
existence and locations of fractal structures. One corol- 
lary from this criterion is that when the initial separation 
velocity is above a certain threshold value, fractal struc- 
tures should disappear — a prediction which agrees with 
our PDE numerics as well as previous numerics on the 
weakly discrete NLS equation (see Fig. 6 in Ref. [27l| ) . 
Lastly, we apply these analytical results to the general- 
ized NLS equations with various nonlinearities such as 
the cubic-quintic, exponential and saturable nonlineari- 
ties, and make detailed predictions on the dynamics of 
their weak interactions. 

This paper is structured as follows. In Sec. 2, we de- 
scribe individual solitary waves in the generalized NLS 
equations. In Sec. 3, we present direct PDE simulation 
results on weak interactions of solitary waves in the gen- 
eralized NLS equations with two different nonlinearities, 
and reveal the common (universal) fractal structures in 
this class of PDEs. In Sec. 4, we analytically derive a 
universal system of dynamical equations (ODEs) for pa- 
rameters of interacting solitary waves using asymptotic 
methods, and show that these ODEs accurately describe 
the weak interactions in the PDEs. In Sec. 5, we solve 
this ODE system analytically when the single parameter 
in this system is equal to zero (which is the integrable 
case). In addition, we derive explicit conditions for the 
solutions of the integrable ODE system to develop finite- 
time singularities. In Sec. 6, we show that fractal struc- 
tures appear in this ODE system when its parameter is 
non-zero, and explore the origin of these fractal struc- 
tures. In Sec. 7, we apply the analytical results to the 
generalized NLS equations with various nonlinearities. In 
Sec. 8, we summarize the results of the paper and make 
some further remarks. 



II. PRELIMINARIES 



The generalized NLS equation is 

iU t + U xx + F{\U\ 2 )U = 0, 



(2.1) 



where F(-) is a real- valued algebraic function with 
F(0) = 0. Equation (|2.1[) supports solitary waves of the 
form 

U = ${x - Vt - Xo -f3)e^ iV ( x - x ^-^ v2t+lf:it - l ' J \ (2.2) 
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where &(6) is a positive function which satisfies the fol- 
lowing equation 



FX|$| 2 )$ - 0$ = 0, 
$ -> 0, |0| -> oo, 



(2.3) 



and (> 0), V, xo, o"o are real constants. For convenience, 
we introduce the notations 



Fi + X0, 

- c , </> 



(2.4) 



Physically, is the propagation constant which is related 
to the solitary-wave amplitude (henceforth, we call an 
amplitude parameter), <fi is the phase of the solitary wave, 
(To is its initial phase, £ is its center position, V is its ve- 
locity, and xq is its initial position. The solitary wave 
is characterized uniquely by its four parameters: V, 0, ctq 
and xq. The asymptotic behavior of this solution at in- 
finity is 



-VP\8\ 



(2.5) 



where c is the tail coefficient which is determined by the 
nonlinear function F and propagation constant 0. Wc 
define the power of the solitary wave as 



P(P) 



$ 2 (6>;0)d<9, 



(2.6) 



which plays an important role in the linear stability of 
the solitary wave. For general functions F, the analytical 
formulas for <!>, P and c are not available. But for some 
special nonlinearities, one can find the analytical solu- 
tions. For instance, for the cubic-quintic nonlinearity 



F(|£/| 2 ) = a|[/| 2 + 7 |£/| 4 , 
the analytical formulas for <f>, P and c are 



450/a 



5 + cosh2V06< : 



P = 



AByfj3{ir/2 - arctan 



n/1 - B 2 



where 



c= y/8B/3/a, 



B = sgn(a)(l 



1607 - 1/2 
3a 2 ' 



(2.7) 

(2.8) 

(2.9) 
(2.10) 

(2.11) 



For special values of a = 1, 7 = 0, Eq. (|2.ip becomes the 
integrable NLS equation, and then 

B = 1, $(0; 0) = v /20sech( v /00), 

P = 4 V / 0, c^v^. (2.12) 



III. UNIVERSAL FRACTAL STRUCTURES IN 
WEAK INTERACTIONS OF SOLITARY- WAVES 

When two solitary waves are placed adjacent to each 
other, they would interfere through tail overlapping. In 
this case, the initial condition is 

U(x,Q) = U 1 (x,0) + U 2 (x,0), 

U k (x,0) = $(x-x o , k ;[3 Q ,k)e i * o - k , 
1, 



t>o,k = -zV Q . k {x - xo k) - fo.fe, 



(3.1) 



where $ satisfies Eq. (|2.3[) . Here "0" in the subscript rep- 
resents the initial value of the underlying parameter. For 
convenience, we assign the left solitary wave with index 
k = 1, and the right solitary wave with index k = 2. 
To study the weak interaction between these two soli- 
tary waves, we require that the two solitary waves are 
both stable, well separated, and having almost the same 
velocities and amplitudes. Introducing notations 

= + 02), V = i(Vi + V 2 ), £ - + &),(3.2) 
and 

A0 = 2 -0 l5 AV = V 2 -V U Ae = 6-fi, (3-3) 

the above requirements then amount to 

Pp > 0, \A(3\ < (3, \AV\ < 1,0A£» 1> | A/3A^| .(3.4) 

Here, Pp = dP/d(3 > corresponds to the Vakhitov- 
Kolokolov criterion for the linear stability of solitary 
waves in Eq. dHJ) [l,[33j]. 

Below, we numerically study the weak interaction of 
solitary waves in Eq. (|2.ip . This equation is numeri- 
cally integrated by the pseudo-spectral method coupled 
with the fourth-order Runge-Kutta integration along the 
time direction. Since each solitary wave has four param- 
eters, we have eight parameters in the initial conditions. 
Due to the phase, translation and Galilean invarianccs 
of Eq. (|2.1[) . we can fix <7o,i = 0, xo.i + xq.2 — and 
Vb = (Vb,i + Vb,2)/2 = without any loss of generality. 
Also, for simplicity, we take AVq = V02 ~ ^0,1 = in all 
our simulations of this section, i.e. the two solitary waves 
are initially at rest. This leaves four free parameters in 
the initial conditions (|3.ip : Axo = xo,2 — xo,i, A<^o = 
0o : 2 - ^0,1 = -00,2, 0o = (/?o,i + A),2)/2, and A/3 = 
00.2 — 0o,i- We define the exit velocity AK» = lim AV. 

t — >+oo 

We also define the collision time t as the time when the 
two solitary waves are the closest (i.e. the separation dis- 
tance between peaks of the two solitary waves the small- 
est) during interactions. The life time of interaction is 
defined as the time length from the beginning (t — 0) to 
the collision time t, which is equal to the collision time 
in value. Thus t will be used to denote the life time as 
well. Of the four parameters in initial conditions, we will 
fix 0o, A0o and Axo, and use A(f>Q as the control param- 
eter and vary it continuously between and 2tt. At each 
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FIG. 1: The graph of exit velocity AVoo versus the initial 
phase difference A</>o in the non-equal initial amplitude case 
of the cubic-quintic NLS equation (12. (|2.7[1 , The cubic and 
quintic nonlinearity coefficients are given in Eq. (|3.5|l . and 
the other (fixed) initial parameters are Axo = 10, /3o = 1, 
A/3 = -0.065, and AVb = 0. 

A0o value, we simulate the evolution of the two solitary 
waves and record the exit velocity and the life time. Nu- 
merically, the exit velocity is determined as follows. We 
let the solitary waves propagate for a long time. If they 
still do not separate, we assign the exit velocity as zero. 
If they do separate, we wait till they have separated far 
apart and their velocities stabilized. Then we locate the 
positions of maximum solitary wave amplitudes at serval 
different time values. The average separation velocity of 
the two solitary waves in these time intervals is assigned 
as the exit velocity. The numerical life time is simply 
the time when the two solitary waves are the closest in 
the simulations. Below, we carry out numerical studies 
of weak interactions as described above on two different 
nonlinearities: the cubic-quintic and exponential nonlin- 
earities. 



A. Weak interactions for the cubic-quintic 
nonlinearity 

Our first example of nonlinearity is the cubic-quintic 
nonlinearity (|2.7[) , which arises in a wide array of physical 
systems such as optics @ and boson condensates [32], [34| . 
In this nonlinearity, we set 

a = 1, 7 = 0.04. (3.5) 

It is easy to verify that all solitary waves (|2.8p in this 
case are linearly stable using the Vakhitov-Kolokolov cri- 
terion. In our simulations, we set Axo = 10 and 0q = 1. 
The x interval is 70 units wide, discretized by 512 grid 
points; and the time step size is 0.004. 

We first study the nonequal-amplitude case, and take 
A/? = —0.065. The AVoo versus A</> diagram is shown 
in FigfJJ The prominent structures in this graph can be 
split into two regions: one region is —0.34 < A0o < 2.5, 



and the other region is 2.9 < A</>o < 3.3. The structures 
in these two regions turn out to be quite similar (except a 
horizontal reflection with respect to a vertical axis), thus 
we focus on the larger region —0.34 < A</>o < 2.5 below. 
The main structure in this region forms a sequence of 
hills; their widths get smaller from the right to the left, 
and their heights are about the same. These hills will be 
called the primary hills. This primary-hill sequence con- 
verges to the accumulation point A<p 0c = —0.339. In or- 
der to see this hill sequence near the accumulation point 
Acj) 0c more clearly, we zoom in the region [—0.35,0.4], 
and the zoomed-in diagram is shown in FigJ5] In this 
figure, the cascading sequence can be seen very clearly 
(see FigJ^a)). In FigOJb), the corresponding life-time 
diagram is displayed. We can see that on the same hill, 
interactions have roughly the same life time. On differ- 
ent hills, life times are different: hills closer to the ac- 
cumulation point A(/)q c have longer life times. Between 
hills, even longer life times can be seen, suggesting more 
complex dynamics there. To explore differences in inter- 
action dynamics on different hills, we select three points 
A0 O = 0.1759,-0.0057,-0.1053 (marked in FigETa) by 
circles) on three adjacent primary hills. These points are 
at the same relative positions (roughly halfway between 
the peak and bottom) of the respective hills. At these 
points, the interaction dynamics is plotted in FigJ^l-3). 
Here only the separation distance A£ versus time t graphs 
are shown. We find that these three dynamical processes 
are similar, except that the oscillation times before final 
separation differ by one from one hill to the next. The 
life times t n of interactions on this primary hill sequence 
are found to be an almost perfect linear function of the 
hill index n as 

ujt n — 2mr + 5, (3-6) 

where the least-square linear fit gives 

lu = 0.08605, S = 2.8897. (3.7) 

Here the life time of each primary hill is measured nu- 
merically at the relative location of that hill shown in 
Figj2ja) by circles. This life-time formula has the same 
form as those for all window sequences reported before 

In addition to the primary hill sequence as described 
above, FigfT] also possesses higher-order structures be- 
tween primary hills. To demonstrate, we first isolate the 
long interval [—0.35, 2.5] in Fig. Q]and re-plot that part 
of the graph in Figj^a). Then we zoom into its sub- 
interval [0.91, 0.995], which is between the two largest 
primary hills in Fig[3ja). The zoomed-in graph is shown 
in FigJ3fb). We see that the zoomed-in graph is simi- 
lar to Fig|3Ja) , but the cascading direction has reversed. 
This behavior is analogous to that reported in [lj| (HI 
for the </> 4 model and the coupled NLS equations. The 
main structure in this zoomed-in window is again two se- 
quences of hills, accumulating to the left and right respec- 
tively. We call them secondary hills. Between secondary 
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FIG. 2: (a) The exit velocity versus initial phase difference graph of Fig. 1 re-plotted near the accumulation point of the 
primary hill sequence; (b) the life time versus initial phase difference graph; (l)-(3): separation versus time diagrams of solitary 
wave interactions at three values of A<fio marked by circles in (a): (1) 0.1759; (2) —0.0057; (3) —0.1053. 
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FIG. 3: Top: the exit velocity versus initial phase difference graph of Fig. 1 and its two zoomed-in structures; bottom: 
soliton-positions versus time diagrams at three values of A</>o marked by circles in the top panel: (1) 1.511; (2) 0.95248; (3) 
0.9669. 



hills, we see even higher-order structures. To see these 
structures more clearly, we zoom into the sub-interval 
[0.9657, 0.96785], which is between the two largest sec- 
ondary hills in Fig[3]Jb). The zoomed-in graph is shown 
in Figj3]^c). We see that it is again similar to Fig. \Mb) 
but with a reversed cascading direction. One can zoom 
into the regions between these tertiary hills in FiglHJc) 
further, and will get even higher order structures which 
are similar to the ones shown in FigEl Thus Fig(3^a) is a 



fractal structure! We have also explored the interaction 
dynamics on this fractal. To demonstrate, we pick three 
A^o values 1.511, 0.95248, 0.9669 which are at the same 
relative positions of the fractal (roughly halfway between 
the peak and bottom of the widest hills) in Figj3{a)-(c) 
(marked by circles). The interaction dynamics at these 
three points are displayed in FigJJHl-3) respectively. Here 
the positions of maximum amplitudes of the interacting 
waves are plotted against time. We see that these dynam- 
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FIG. 4: The exit velocity versus initial phase difference 
graph in the equal initial amplitude case of the cubic-quintic 
NLS equation 12. ip . (|2,7[l . The cubic and quintic nonlinearity 
coefficients as well as the initial conditions are the same as in 
Fig. f , except that Aflo = now. 

ical patterns are clearly similar, except that the numbers 
of oscillations before final separation are different. 

In the above numerical simulations, the two solitary 
waves have different initial amplitudes (A/3o = —0.065). 
We have also studied interactions of equal-amplitude soli- 
tary waves, i.e. with A/3o = 0, while keeping the other 
parameters the same. In this case, the graph of exit ve- 
locity AVqo versus initial phase difference A</> is shown 
in Figf?] This graph is symmetric with respect to A</> for 
obvious reasons. Examination of this graph shows that 
it is also a fractal. Thus fractal dependence arises in 
weak interactions of both equal and non-equal amplitude 
solitary waves. 



3n/2 



0.; 



AV o. 




FIG. 5: The exit velocity versus initial phase difference graphs 
for the exponential nonlinearity 1)3.80 : (a) the non-equal ini- 
tial amplitude case with A/3q = —0.045; (b) the equal initial 
amplitude case with A/3o = 0. The other (fixed) initial pa- 
rameters are /3o = 2.3, Aa;o = 8, and AVb = 0. 



B. Weak interactions with exponential nonlinearity 

To explore whether the above fractal structures for 
weak interactions persist or not with other types of non- 
linearities, we consider in this subsection a different type 
of nonlinearity — the exponential nonlinearity, with 

F(\U\ 2 ) = e m - 1. (3.8) 

Here, —1 is introduced into this function to meet the 
condition F(Q) = 0. Note that this nonlinearity does not 
have any parameters. Throughout this subsection, we set 
the initial separation Axq = 8, and average propagation 
constant /3q = 2.3. We study two cases, one for non-equal 
amplitudes with A/?o = —0.045, and the other for equal 
amplitudes with A/?o = 0. For both cases, the control 
parameter is A0o as before. In our simulations, the x in- 
terval was 70 units wide, discretized by 512 grid points. 
The time step size was 0.002. The AV^ verse A(f>o graphs 
for both cases are plotted in Fig. \5\ We have verified that 
both graphs in this figure are fractals. Comparing these 



fractals with those in Figs. Q] and |4] of the cubic-quintic 
nonlinearity, we see that the fractal structures for these 
two different nonlinearities are very similar. The only 
major difference between them is in the non-equal ampli- 
tude case, where there is only one primary hill sequence 
(accumulating toward the left) for the exponential non- 
linearity, while there are two primary hill sequences for 
the cubic-quintic nonlinearity. It is remarkable that two 
very different nonlinearities exhibit quite similar fractal 
dependence on initial conditions. Thus fractal scattering 
appears to be a universal feature in weak interactions of 
Eq. 1)2. ip rather than an accident. This leads us to the 
following questions: how can we analytically establish the 
universal nature of fractal scatterings for Eq. (|2.1[) with 
general nonlinearities? how can we analytically explain 
the major differences of fractals for different nonlineari- 
ties? These questions will be answered in the following 
sections. 



7 



IV. DYNAMICAL EQUATIONS 

To study weak interactions analytically, we use the 
Karpman-Solov'ev method [§] by treating the interfer- 
ence as a small perturbation to each solitary wave (see 
also (HI). This method has been successfully used be- 
fore on the NLS equation @, H, Ell the modified 
NLS equation [llj, the Manakov equations [13], as well 
as the (non-integrable) coupled NLS equations [36[. To 
proceed, we first need to consider the evolution of a single 
solitary wave in the perturbed generalized NLS equation 



iU t + U xx + F(\U\ 2 )U = eG, 



(4.1) 



where function G is a perturbation term, and e is a small 
parameter. Without perturbations (e = 0), the solitary 
wave (|2.2p is an exact solution of Eq. (|4.ip , and its internal 
parameters V,j3, aoi x o are time-independent. When the 
perturbation is turned on, these internal parameters of 
the solitary wave will evolve slowly on the time scale T = 
et. The multiple-scale perturbation theory for this slow 
evolution is well known [3(| [37j • We write the perturbed 
solution as 



where 



U = <f>(9,t,T)e 



Vdt — xq , 



iVB/2+ia 



(4.2) 



(f3 + V 2 /4)dt-a a . (4.3) 



Here V(T), /3(T), a (T), x Q (T) are all functions of slow 
time T. Next, we will derive the dynamical equations 
(ODEs) for the slow-time evolution of these parameters. 
Substituting (14. 2|) into (|4.1j) . we get the equation for $ 

as 

i$t + l> ee -/3l> + F($ 2 )$ = 



-e (Vx 0T /2 - V T e/2 + a 0T ) 



(4.4) 

where <j> is defined in Eq. (|2.4p . We expand the amplitude 
function $ into a perturbation series 



$ = $(0; (3) + e$ + 0(e 2 



(4.5) 



The equation at order e° is satisfied automatically since 
$ satisfies Eq. (|2.3j) . At order e, the equation for $ can 
be written as 



where 



* = 



Co 



W t = H, 



$ + <r \ _ / o £ 
^ I A o 



-do 
-d„ 



+ 13- 
+ 0- 



F($ 2 ) 
_F($ 2 ) 



2$ 2 f ($ 2 



(4.6) 



(4.7) 



(4.8) 



and 
H = 



-G*e^ + Ge-** - 2i§0 T + 2i<S> g x 0T 

- Ge-^ + (Vxq T - 9V T + 2a 0T )^ 



(4.9) 



Here the superscript " *" represents complex conjugation. 
Operator C has two eigenfunctions and two generalized 
eigenfunctions associated with the zero eigenvalue, 



*i = 
*i = 



$0 




■6$/2 



*2 = 



*2 = 




4> 





with the relations 
C^! k = 0, 



(4.10) 



(4.11) 



In order for the inhomogeneous solution '5 of the first- 
order equation (|4.6[) to be non-secular at large time, the 
inhomogeneous term in Eq. (|4.6|) must be orthogonal to 
the above eigenfunctions and generalized eigenfunctions 
of the zero eigenvalue, i.e., 



{H, * fc ) = {H, V k ) = 0, k 
under the inner product defined as 



1,2, 



dO. 



(4.12) 



(4.13) 



Here is the Hermitian of Ff.. Evaluating the four 
integrals in Eq. (|4.12j ). the slow-time evolution equations 
for parameters V(T), /3(T),a (T),x (T) will be obtained. 
These evolution equations can be written as 



P— = 2 / SflfGV* + Ge~**)d0, 
dT J_ 0o 



(4.14) 



e lrt >)d6, (4.15) 



P 



dxo 
~dT 



$6>(Ge 



-ic£> 



G*e**)d0, (4.16) 



, V dxo 



j /■oo 

Pp{1 2 W + ^F ) = J 00 + Ge-»)M{A.17) 

These equations will be critical for the development of 
weak interaction theory of solitary waves below. 

Now, we consider the weak interaction of two solitary 
waves. Here the tail overlapping can be considered as a 
small perturbation which causes the internal parameters 
of each solitary wave to evolve on a slow time scale et. 
Here e is the magnitude of tail overlapping which is ex- 
ponentially small with solitary wave spacing A£. We will 
not introduce e explicitly in the next analysis. To the 
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leading order, the interacting solution is simply a super- 
position of two solitary waves, 

U = Ui + U 2 , 
U k = fc = l,2, (4.18) 

where all parameters slowly vary over time. Picking 
up the dominant interference terms, each solitary wave 
is governed by the following perturbed generalized NLS 
equations: 

%U k , t + iU kM + F{\U k \ 2 )U k = H k , (4.19) 

where 

H k = - (F(\U k \ 2 ) + F'(\U k \ 2 )\U k \ 2 )U^ k 

- F'(\U k \ 2 )U 2 U^ k . (4.20) 



In this paper, we only study the weak interaction, so 
conditions (|3.4p are assumed. Since |AV| <C 1, the phase 
difference 

A(f> = fa - <f>i re -VA£/2 + Act, (4.21) 

which is independent of 9. 

Now We apply the above solitary wave perturbation 
theory to Eq. (|4.19[) . In this problem, 

eGe-<+= - (F($l)$^ k + F\$l)$l<t> 3 _ k )e(-V k+liA * 
- F , ($l)Ql*a- k e<-- r > hi *+. (4.22) 

Substituting (|4~22"1) into Eqs. (|4~14)l - ([4~T7| . we obtain the 
following dynamical equations 



Pk^- = -iJ $ M (F($ 2 )$ 3 -fc + 2F'($l)$^ 3 _ fc )^cos(A^), 
Pk,h^ = (-i) fe 2 £ $ h F(*\)$ i - h d0MN>), 



(4.23) 



(4.24) 



dx ky o 
1 dt 



P k,/3 k ( 



V k dx k n , da k n . 



2 dt 



dt 



/CO 
®k0F($l)$ 3 ^ k d8 sin(A^), 
-00 

/OO 
$ kA (F(^)4 3 -fc + 2F'($l)$l$ 3 - k )d6 cos(A0), 
-OO 



(4.25) 



(4.26) 



where P k , k = 1,2 are powers of the two individual soli- 
tary waves. These equations can be simplified greatly. 
Due to assumptions (|3.4p . and noticing that $(#) and 
$p(0) are even functions of 6, the leading-order terms 
of the above integrals can be explicitly calculated. For 
instance, 



Similarly, 

/OO 
-oo 

/CO 
$ e (F($ 2 ) + 2F'($ 2 )$ 2 )ce(- 1 ) fc+1 ^ 8 ^e-^ JA « 
-OO 

/CO 
$F($ 2 )ce^ e ^e-^ A « 
-oo 



(-l) fe 2/3c 2 e 



,2 -V?A? 



(4.28) 



$ fc F($ 2 )$ 3 -fc^ 
f° <S>F{$> 2 )ce^ k+1 ^ e dee~^ 

J — OO 

/°° (/?$ - $ ee ) Ce ^^e-^ A « 

i/ oc 



(4.27) 



$ fe 0F($ 2 )$ 3 -fc^ 

/OO _ _ 

<S>9F(<S> 2 )ce^ k+1 ^ e dee-^ A Z 
-oo 



oc 

OO 



(-l) fc 



fc+1 



$0F{$ 2 )ce^ e d0t 



(4.29) 
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/CO 
-co 

/CO 
^(JX* 2 ) + 2F'($ 2 )$ 2 )ce(- 1 ) fc+1 ^^e-^ A « 
-co 

/CO _ 
$,g(F($ 2 ) + 2F'($ 2 )$ 2 )ce v?e d6le- v?A « 
-co 

= f D a e~^. (4.30) 

With the above simplifications, the dynamical equations 
reduce to 

P^t = (-i)fe+i 8 /3c 2 cos(A^)e-^ A «, (4.31) 
at 



P^ = (-l) fe 4 v ^c 2 sin(A0)e-^ A «, (4.32) 



p dxkfl = _ 2Di sin(A0)e -v^A ?j (4 33) 

at 



P ^V~^r + ^df ) = ~ 2D * ^(A^e-^, (4.34) 

where P is the power of the solitary wave with propaga- 
tion constant ft and D\,D 2 are defined in (|4.29[) . (|4.30|) . 
From the above equations, we find that 



ft = V t = 0, 



A& = Ay, 



= A/3, 



Aft = 



!6ft 2 



8Vft 2 



■ cos(A</>)e~ 



sin(A0)e-^ A5 . 



(4.35) 



(4.36) 



(4.37) 



(4.38) 



(4.39) 



Equations (|4.35[) - (|4.39[) are the key results in the weak 
interaction theory of solitary waves. These equations can 
be further simplified by variable rcscalings. Introducing 
notations 



V = A0,C = -v^Ae,/ 



16ft/ 



2J2 



Pa 



,(4.40) 



and 



P 



then the dynamical equations (|4.36p - (|4.39fl reduce to 



(4.41) 



{ 



Q TT = cos V>e 

iPtt = (1 + e) sin-0e < » 



Eq. (|4.42p is the final dynamical system we obtained 
for the analytical treatment of weak interactions in the 
generalized NLS equations ()2. ljl . It is important to re- 
mark that Eq. (|4.42[) is universal for the generalized NLS 
equations with arbitrary nonlinearities. It contains only 
a single parameter £, which depends on the specific form 
of nonlinearity. 

Eq. (|4.42|) has the following general properties. First, 
it is Hamiltonian with the conserved Hamiltonian (en- 
ergy) as 

£ = ^(C 2 -^ 2 )-e c cosV+ * J 2 - (4.43) 
2 2(l + e) 

Here ( ) = d/dr. Second, it has some symmetry proper- 
ties. One is that it is time-reversible, i.e., if [C(t), V , ( r )] 
is a solution with initial conditions (Co, Co, V'o, V'o), then 
[C( — t), V( — t)] is a solution with initial conditions 
(Co, — Co, "007 — V'o)- Another symmetry is on phase flip- 
ping, i.e., if [C(t),iP(t)] is a solution with initial condi- 
tions (Co, Co, V'o, V'o), then [C(t), — ip( T )) is a solution with 
initial conditions (Co, Co, "V'o, "V'o)- Physically, this lat- 
ter symmetry corresponds to the interchange of the left 
and right solitary waves in the PDE evolutions, which 
of course does not change the interaction outcome. Eq. 
(14.42[) also has the property that if is a solution, so 
is iP(t) + 2mr for any integer of n. This reflects the fact 
that in the PDE system, solution evolution remains the 
same if the phase difference between the solitary waves 
changes by a multiple of 2ir. 

The dynamical equations (|4.42p are asymptotically ac- 
curate in describing weak interactions in the PDE system 
(to the leading order) when the spacing AC is large. Sur- 
prisingly, even when the two solitary waves come close to 
each other, Eq. (|4.42[) can still describe the interaction 
process very well. This is analogous to weak interac- 
tions in the (integrable) NLS equation Below, we 
make detailed comparisons between the ODE solutions 
of Ea. p~32"f and the PDE solutions in Sec. 3(A) for the 
cubic-quintic nonlinearity. Inserting the parameter val- 
ues a = 1,7 = 0.04 and ft = 1 of PDE simulations into 
Eqs. (EU) and (j2"3U]) . we get P = 3.74720, Pp = 1.64835, 
and c = 2.69495, thus / = 31.01080, g = 35.24845, and 
e = 0.13665. Corresponding to the initial conditions for 
PDE simulations in Sec. 3(A), the initial conditions for 
the ODE system (I4.42[) with nonequal and equal initial 
amplitudes are 



Co = -10, Co = 0, V'o = -0.01167, 



and 



Co 



-10, Co = 0, Vo = o, 



(4.44) 



(4.45) 



(4.42) 



respectively. In both cases, the initial phase difference 
V'o is the control parameter as in PDE simulations. The 
ODE system (|4.42l) is numerically solved by the fourth- 
order Runge-Kutta method, with the time step set as 
0.01. The simulation results on the exit velocity —Coo 



10 



0.14 



0.14r 




FIG. 6: The exit velocity (— Coo) versus the initial phase difference (ipo) graphs from the ODE model (14.42[) . Here the initial 
conditions for (a) and (b) are chosen corresponding to the PDE simulation results in Figs. 1 and 4 respectively (see text). 



versus tpo graph are shown in Fig[5J Clearly, these graphs 
are very similar to Figs. [T]and|4]) from PDE simulations. 
We have also investigated the detailed structures of Fig. 
[6] in ways analogous to what we did for Figs. 2 and 3. 
Specifically, we have examined the primary hill sequence 
in Fig. E^a), and zoomed into regions between primary 
hills. The results are shown in FigsH and [8] respectively. 
Both figures closely resemble Figs. 2 and 3 from the PDE 
simulations. 

The agreement between the ODE model and the PDE 
simulations is not only qualitative, but also quantitative. 
To demonstrate, we compare the locations and life times 
of primary hill sequences in Figs. 2 and [7] The com- 
parison results are summarized in Table 1. Very good 
quantitative agreement between them can be seen. In the 
ODE model, the life time is also an almost perfect linear 
function of the hill index n in the form (|3.6p . When the 
time rescaling (|4.41l) is recovered, the ODE model gives 

u\ode = 0.08570, S\ode = 2.9655, (4.46) 

closely resembling the corresponding values (|3.7p from 
the PDE simulations. 

Above we have established that the reduced ODE sys- 
tem (|4.42p accurately describes weak interactions of the 
PDE system. Since the ODE system (I4.42[) is universal 
for Eq. (|2.1[) regardless of details of its nonlinearities, 
we see that the hill sequences and fractal structures in 
Eq. (|4.42p are universal for weak interactions of solitary 
waves in the PDE system (|2~Tj) . as Figs. HJ IH andll 
clearly indicate. 

Next we will turn our attention to the ODE system 
(|4.42p . and analyze its solution dynamics in more detail. 
In particular, we would like to understand why fractal 
structures arise in this system, and how to analytically 
predict their locations and other main features. 



TABLE I: Comparison on locations and life times of primary 
hills in Figs. 2(a) and 7(a) from the PDE and ODE simula- 
tions. 



location (PDE) location (ODE) life (PDE) life (ODE) 



n 


A<f>0,n 


V>0,n 


t n 




1 


1.7735 


1.7794 


65 


68 


2 


0.6985 


0.7015 


117 


119 


3 


0.2430 


0.2468 


183 


185 


4 


0.0280 


0.0359 


253 


255 


5 


-0.0850 


-0.0763 


325 


327 


6 


-0.1530 


-0.1431 


398 


400 


7 


-0.1963 


-0.1863 


470 


473 


8 


-0.2258 


-0.2157 


544 


547 


9 


-0.2475 


-0.2367 


617 


620 


10 


-0.2630 


-0.2523 


691 


693 


CO 


-0.3392 


-0.3280 


oo 


oo 



V. SOLUTIONS OF THE INTEGRABLE 
DYNAMICAL EQUATIONS AND THEIR 
SINGULARITY CONDITIONS 

Eq. (l4.42p conserves energy (|4.43|) for all values of e. 
When e = 0, it has another conserved quantity, 

M = CV> ~ e c sinV', (5.1) 

which can be called the momentum of this system. In this 
particular case, system (|4.42l) is an integrable Hamilto- 
nian system and can be analytically solved. Let Y = 
C + iip, Eq. (|4.42p becomes 

Y TT = e Y . (5.2) 

The general solution of this equation is 

Y(t) = In [-2C 2 sech 2 (Cit + C 2 )] , (5.3) 
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FIG. 7: (a) The exit velocity versus initial phase difference graph of Fig. [6ja) re-plotted near the accumulation point of the 
primary hill sequence; (b) the life time versus initial phase difference graph; (l)-(3): separation (— C) versus time (r) diagrams 
at three values of tpo marked by circles in (a): (1) 0.187; (2) 0.0056; (3) —0.0939. All these graphs are obtained from the ODE 
model (|4.42fl . and they should be compared to the corresponding PDE graphs in Fig. [2] 




XXX 



FIG. 8: Top: the exit velocity versus initial phase difference graph of Fig. [6] and its two zoomed-in structures; bottom: 
separation versus time diagrams at three values of i/>o marked by circles in the top panel: (1) 1.532; (2) 0.95071; (3) 0.96603. 
These graphs from the ODE model should be compared to the corresponding PDE graphs in Fig. [3] 



where 



and 



C x = \ v/y 2 -2e^ = -L VETiM, (5.4) 



c 2 



-arctanh 



Y 
2Ci 



(5.5) 



Y which differ by a multiple of 2tti correspond to the 
same physical solution, thus it does not matter which 
Riemann surface one takes for the logarithmic function 
in Eq. (JOJ). If d = 0, i.e., Y = ±V2e Yo/2 , the solution 
(15. 3p degenerates to the form 



Here the branch of the square root function in (|5.4[) is 
chosen such that Re(Ci) > 0. It is noted that solutions 



Y(t) = -2 In 



,-hYo 



(5.6) 
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The asymptotic behaviors of these solutions as r — > oo 
can be easily determined. Let 



C\ = a + bi, 



CI 



c + di, 



(5.7) 



where a, 6, c, d are real constants, then the following 
leading-order asymptotic expressions for the solution can 
be obtained when r — > oo: 

(1) a^O: Y(t) -2|a|r-sgn(a)26ri; (5.8) 

(2) a=0, 6^0: 
(2a) d = : 

F(r) = In 26 2 - In cos 2 6(r + c); (5.9) 
(26) d ^ : 

Y(t) = In 46 2 - In [cosh 26d + cos 26(r + c)] 

+2i arctan [tanh(fcd) tan6(r + c)] ; (5.10) 

(3) a = 0,6 = 0: F(r) -> -21n(=Fr). (5.11) 

From these asymptotic expressions, we see that when 
a 7^ 0, the two solitary waves eventually move away 
from each other with exit velocity 2\a\; when a = but 
6^0, the solution is time-periodic for both d = and 
d ^ 0, the difference being that in the former case, the 
periodic solution exhibits finite-time singularities (where 
C = Re(Y) — * oo), while in the latter case, the solu- 
tion has no singularities; when a = b = 0, the two soli- 
tary waves eventually separate logarithmically, and the 
exit velocity is zero. As an example, we take the ini- 
tial conditions (|4.44|) . In this case, the graph of exit 
velocity — Coo(= 2|a|) versus V'o is plotted in Fig. [5] (bot- 
tom panel). This graph is smooth everywhere, except 
at ipQ = 0, ±7r where it has a cusp (due to the absolute- 
value function in \a\). The squares and diamonds on this 
graph will be explained later. Clearly, this graph has no 
fractal structure anywhere. Thus, fractal dependence is a 
signature of the dynamical system (|4.42|) when it is non- 
integrable (with e ^ 0), not when it is integrable (with 
e = 0). 

The above asymptotic states do not tell the full story 
about the solution dynamics in the integrable system. 
For instance, for the case of a ^ 0, even though the so- 
lution has a benign- looking asymptotics (|5.8[) as r — ► oo, 
the solution can still develop a singularity (where the sep- 
aration £ — > oo) at a finite time. These solutions with 
finite-time singularities turn out to be critical for the ap- 
pearance of fractal structures in the non-integrable sys- 
tem, as our numerics in the next section will indicate. 
Thus we analyze these singularity solutions in more de- 
tails below. The necessary and sufficient conditions for 
singularities in solution (|5.3|) are that 



(5.12) 



cosh(Cif + C 2 ) = 0, 



and f > 0, where f is the time of singularity. If f < 0, 
i.e., singularities in the solution occur at a negative time, 
such singularities are irrelevant for the time evolution of 



Eq. (|4.42[) and need not be considered. The solutions of 
Eq. (HTTP are 

Cxf + C 2 = -(2n + l)m, n = 0, ±1, ±2, • • • . (5.13) 

This is a complex-valued relation, which gives two real 
relations on f, C\ and C 2 . When a ^ 0, i.e., C% is 
not purely imaginary, we find by separating the real and 
imaginary parts of Eq. (|5. 13|) that the solution (|5.3[) has 
a single finite-time singularity of the type ln(r — f) if the 
following conditions are satisfied: 



S 



Im(CTC2 
Re(Ci) 



1 



(2n+l)7r, 



0,±1,±2, 



Re(C 2 
Rc(Ci 



>0. 



(5.14) 



(5.15) 



Here Re(-) and Im(-) represent the real and imaginary 
parts of a complex number. When a — (b ^ 0), sin- 
gularity solutions exist if d = 0. These solutions have 
an infinite number of finite-time singularities of the type 
ln(r — f), as the formula (|5.9|) indicates. Physically, at 
the time of singularity f, the two solitary waves strongly 
collide, thus f is the collision time. Whether conditions 
(15.14[) and (|5. 15[> can be satisfied depends on the ini- 
tial conditions (which determine the C\ and C 2 values, 
see (|5.4I) . (I5.5[) ). In the text below, we will call initial 
conditions (Co, Co; V'O; V'o) which satisfy Eqs. (|5.14|) and 
(|5.15|) as singularity points. At singularity points, solu- 
tions of the integrable dynamical system (|4.42|) develop 
finite-time singularities. 

To demonstrate how to determine singularity points 
in the initial-condition space, we take initial conditions 
(I4.44p of Fig. HJa) as an example. Here i[>o is a control pa- 
rameter. With these initial conditions, the graph of func- 
tion S(ipo) is plotted in Fig. [9] (top panel). This graph 
has a maximum 0.96. As tpo — > + or 7r~, S(ipo) — > — oo. 
As we can see from this graph, for any value of n < — 1, 

Eq. (|5.14p has two roots, ^q 1 ^ and ip(^n- We have checked 
that these roots satisfy the other singularity condition 
(|5.15[) . thus these an d values are singularity 
points. It is noted that the graph of function S(ipo) also 
has another piece in the interval n < tJjq < 27r, which 
is the mirror image of that shown in Fig. [9] around the 
point "00 = t- But in that interval, f < 0, not satisfying 
the second singularity condition (|5.15[) . thus we did not 
plot that piece of the graph in Fig. [9l 

Next, we examine these singularity points and 

(2) 

ipQ n in more detail. These points form two infinite 
sequences with n = —1,-2,..., which accumulate at 
ipa = + and 7T~ respectively. In Fig. [5] (bottom panel), 
these two sequences are marked by squares and diamonds 
on the exit velocity versus ipo graph. Calculating the 
asymptotics of C\ from Eq. (I5.14|) and substituting it 
into Eq. (|5. 15|) , we find that the collision times f n of both 
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FIG. 9: Top: graph of the function S(tpo) defined in Eq. 
(|5.14p for the initial conditions (I4.44|l . Intersections of the 
graph with horizontal lines are singularity points. Bottom: 
exit velocity versus V>o graph in the integrable system (|4.42l l. 
Both squares and diamonds are singularity points. 



sequences have the following asymptotic expressions: 

ujf„ — 2\n\ir + 7r, n — > — oo, (5.16) 

where ix> = 2Im(Ci|^, 0= o) and 2Im(Ci |^, 0=7r ) for the left 
and right sequences respectively. The asymptotic formu- 
las for the locations of these two singularity sequences 

{^Onl an d {V'on} can a l so ^ e calculated. We find from 
(|5TT3|| that 



and 



where 



and 



/(2) 



(1) 



At 



(2n+ 1)tt' 



A 2 



(2n+ 1)tt 



?1 — * — OO, 



n — > — oo, 



(5.17) 



(5.18) 



A t = 8e^ Co Re(C 2 )Im 2 (Ci; 



A 2 = 8e- Co Re(C 2 )Im 2 (Ci; 



The above detailed analysis on singularity points was 
performed for the particular initial conditions (|4.44D 



where the two solitary waves are initially stationary 
(Co = 0). What will happen if Co 0? To answer this 
question, we fix Co and ipo as in Eq. (|4.44[) . vary Co, 
and examine how singularity points move in the (tpo, Co) 
plane. The results are shown in FigJTOJ The top curve 
corresponds to n = —1, the next curve corresponding 
to n — —2, and so on. All curves are bounded from 
both above and below except the top one (with n = — 1). 
When n — > — oo, these curves approach the accumula- 
tion curve plotted by the dashed line in FigOjl] Below 
this accumulation curve, there are no singularity points. 
The analytical formula for this accumulation curve can 
be easily derived. On this accumulation curve, C\ must 
be pure imaginary, thus 



M = Coc^o 



sin^o = 0, 
e Co cos < 0. 



(5.19) 
(5.20) 



Here (Coc V'o) is an accumulation point. From Eq. (|5.19p . 
we see that the function of the accumulation curve is 



sin -00 



0c 







The maximum and minimum of this curve are 



Coc 



"Coc 



(5.21) 



(5.22) 



For the Co and ip values of Fig. [Hjl we get Coc.min = 
—0.00389. If Co < Coc.minj there are no singularity solu- 
tions for any value of ip . When Coc.min < Co < Coc.max, 
two infinite sequences of singularity points can be found. 
When Co > Coc.max, however, the number of singularity 
points becomes finite; this number gradually decreases 
(down to one) as Co increases. 

The above calculations of singularity points and their 
accumulation curves were made for special choices of ini- 
tial conditions Co = —10 and tpo = —0.01167 (see Fig. 
[TO]). In view of the importance of singularity points for 
fractal structures which we will reveal in the next sec- 
tion, we would like to discuss singularity points and their 
accumulation curves further for general initial conditions 
below. 

First, we examine the accumulation curve in the (^o, 
Co) plane for general initial conditions Co and ipQ. In this 
general case, the accumulation curve (if it exists) is nec- 
essarily given by Eq. (|5.21[) . But the curve (|5.21|) (or por- 
tions of it) may not satisfy condition (|5.20[) . thus may not 
actually be the accumulation curve. Below we determine 
what portions of the curve l|5.2ip are the accumulation 
curve. Before we do so, let us first point out that con- 
ditions (|5.19[) and (|5.20p are not only the necessary, but 
also sufficient conditions for the accumulation curve. In 
addition, the accumulation of singularity points toward 
the accumulation curve is always from the upper side, 
not lower side. To show these, we only need to prove that 
condition (|5.15p holds only on the upper edge of the curve 
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FIG. 10: Singularity points satisfying conditions (I5.14|l - (|5.15|) 
in the (0o, Co) plane. The dashed curve is the accumulation 
curve. Here Co = -10, O = -0.01167. 



(|5.2ip . but not the lower edge of it. On the upper edge 
of (|5.21[) . C\ is purely imaginary, and sgn(M) = sgn(^o)- 
Thus sgn(Im(Ci)) = sgn(M) = sgn(0 o ). Consequently, 
Re(Yo/2Ci) > 0. Notice that for any complex num- 
ber z, Re[tanh(z)] and Re(z) have the same sign, hence 
Re(C 2 ) < 0. Then due to Re(Ci) > 0, condition (|5~l"5|) 
thus holds. By similar reasoning, we can show that on the 
lower edge of the curve (|5.21[) . condition (|5.15[) does not 
hold. Thus singularity points accumulate toward (|5.2ip 
only from above, not below. 

Now we turn to equations (|5.19p and (|5.20p , and use 
them to determine the accumulation curve for the general 
case. Substituting Eq. (|5.2ip into inequality (|5.20p and 
simplifying, this inequality becomes 



^ ($j + e c °(l + cosV-o)) (i>l - e Co (l - cos0 o )) < 0, 

(5.23) 



which is equivalent to 

cosset > 1 



(5.24) 



Thus the accumulation curve is the parts of curve (|5.21|) 
where 0o satisfies the constraint (I5.24[) . If ipo > 2e < =°, 
condition (|5.24[) is satisfied for all values of i/jq, hence the 
entire curve (|5.2ip is the accumulation curve. If < ipo < 
2e^° , portions of the curve (|5.2ip centered at ipo = tt do 
not satisfy condition (|5.24p . thus do not belong to the 
accumulation curve. The rest of the curve (|5.2ip does 
satisfy condition (|5.24p . thus is the accumulation curve. 
If ipo — (equal initial amplitude case), no value of 0o 
satisfies condition (|5.24p . thus accumulation points do 
not exist. 

Next, we derive two general properties about singu- 
larity points in the (t/'o, Co) plane for general initial con- 
ditions Co an d 0o- One property is that, if Coc is on 



the accumulation curve, then for any Co < Coc singular- 
ity points can not exist. We will prove this by showing 
that f < for Co < Coc- To show f < 0, we only need 
to show Re(Y /2Ci) < (see above). Without loss of 
generality, we only show this for ipo < 0; the proof for 
"00 > is similar (in fact, as has been pointed out be- 
fore, flipping the sign of ipa physically amounts to inter- 
changing the positions of the left and right solitary waves 
and thus does not affect the interaction outcome). For 
ipo < and Co < Coc Ci is in the first quadrant (as 
M > 0). If Co < 0, then Y is in the third quadrant, thus 
Re(Y /2Ci) < holds. Now we consider < Co < Coc- 
In this case, Yq is in the fourth quadrant, hence iYo lies 
in the first quadrant (like Ci). To show Re(Y /2Ci) < 0, 
we only need to show arg(iY ) < arg(2Ci). Since both 
iYq and C\ are in the first quadrant, we only need to 
show arg(-Yp) < arg(4C 2 ). Notice that 



- Y l 



4C 2 = -2e Yo , 



(5.25) 



which is independent of Co- In addition, the angle of 
-2e Y ° falls in between those of -Yq 2 and AC' 2 . Thus 
to show arg(— Y 2 ) < arg(4C 2 ), we only need to show 
arg(-Y 2 ) < arg(-2e Yo ). Note that 



Yo = $ - Co 2 



2#oCo, 



(5.26) 



whose angle is an increasing function of Co when 0o < 0, 
thus for Co < Coc 



arg(-Y 2 ) < arg(^ 



t 2 

SOc 



2iV»oCoc)- 



(5.27) 



Now recall that Coc lies on the accumulation curve, thus 
it satisfies the conditions (|5.19[) and (|5.20p . Substituting 
these conditions into (|5.27p . and recalling our assump- 
tions of ipo < and < Co < Coc we find that the right 
hand side of (|5.27p is less than arg(— 2e Y °), thus inequal- 
ity arg(— Yg) < arg(— 2e Y °) is proved. Summarizing the 
above arguments, we conclude that for any Co below the 
accumulation curve, singularity points do not exist in the 
(0o, Co) plane. 

Another general property about singularity points is 
that, at sufficiently large values of Co, there is a unique 
singularity point in the tpo interval. The proof is as fol- 
lows. It is easy to check that when Co + Wo 3> and 
I Co I ^ I Wo I, functions S and f have the following leading- 
order asymptotic expressions, 



l ln 2(C o 2 +0 o 2 ) 
Co eCo 



s 2 S gn(Co)V'o + S A , 



(5.28) 



(5.29) 



where 



0o , 2(C o 2 +0g) 0o 1 ,, qn . 

Sa = — — m ^— — - arctan— n, 5.30 

2|Co| e& | Co | 2 
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and the relative errors are O ^e^° /((g + ipo)j ■ From Eq. 

(|5.29p . we see that the rise of S value over the interval 
< ipo < 27r is 7r, which guarantees that Eq. (|5.14[) has 
a single solution in the tpo interval for a single value of n. 
From Eq. (|5.28[1 . we see that when Co > is sufficiently 
large, f > over the entire tp interval. Thus singularity 
conditions (|5.14[) and (|5.15p admit a unique singularity 
point. We note by passing that when Co is sufficiently 
large negative, f < over the entire tpo interval, thus 
there can not be any singularity points. This is consistent 
with the previous general property proved above. 

To summarize the above results on singularity points 
and accumulation points and slightly extend them, we 
present the following classifications on singularity solu- 
tions in the integrable system (|4.42[) : 

1. tpo = (equal initial amplitudes): 

• If Co > — v / 2e'» / 2 , a singularity solution exists 
at the single singularity point tpo = 0. Here 
M = 0, and E > for Co > \/2e^°/ 2 and E < 
otherwise; 

• If Co < -V2V / 2 , there are no singularity so- 
lutions for any ipo', 

2. tpo 7^ 0, Co = (non-equal initial amplitudes, zero 
initial velocities): 

• If ipQ > 2e'> , singularity solutions exist at two 
infinite sequences of tpo values, accumulating 
at -00 = {0+,7r~}, or {7r+,27r~}, for tpo < 
and tpo > respectively; 

• If < ipQ < 2e^° : singularity solutions exist at 
one infinite sequence of tpo values, accumulat- 
ing at tpo = + or 2ir~ for tpo < and tpo > 
respectively; 

On these sequences of singularity points, M 7^ 
and E =/= generically (at the accumulation points, 
M = 0, and E < 0); 

3. tpo 7^ (the general non-equal initial amplitude 

case): 

In this case, the accumulation curve is the parts of 
curve (|5.2ip where tpo satisfies the constraint (|5.24|) . 
When tpQ > 2e < "°, the entire curve (|5.21[) is the ac- 
cumulation curve. When <ipo < 2eCo > the 

accu- 
mulation curve is l|5.2ip except portions of it which 
are centered at tp = 7r. 

For Co below the accumulation curve, there are no 
singularity points; at sufficiently large Co values, 
there is a single singularity point. 

At all these singularity values, E and M are non- 
zero generically (except the accumulation points 
where M = 0). 
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FIG. 11: The exit velocity versus initial phase difference 
graphs in the ODE model (|4.42|l at various values of e: (1) 
0.13665; (2) 0.036; (3) 0.0036; (4) 0; (5) - 0.0036; (6) -0.036. 
The initial conditions are given in (I4.44|l . The squares and 
diamonds in (4) are singularity points of the integrable system 
(see Fig. 9, bottom). 



It is noted that in the above classifications, case (2) is 
just a special case of case (3), and can be readily deduced 
from (3). Case (1) can be deduced from (3) as well under 
the limit tpo — > 0. But cases (1) and (2) are important 
special cases, hence we listed them out separately. 



VI. ORIGINS OF FRACTAL STRUCTURES IN 
THE NON-INTEGRABLE DYNAMICAL 
EQUATIONS 

We have known from Figs. O [7] and [5] that the non- 
integrable ODE system (|4.42[) exhibits hill sequences and 
fractal structures which coincide with those in the PDE 
simulations, but such structures do not exist when this 
ODE system becomes integrable. The natural question 
then is: where do the fractal structures in the non- 
integrable system (|4.42l) come from? In this section, we 
will establish through careful numerics that these fractal 
structures bifurcate from singularity points of the inte- 
grable system. 

To determine the origin of these fractals, we take the 
same initial conditions (|4.44|) as in Fig. [|3a), but grad- 
ually decrease the value of e from 0.13665 of Fig. [Ha) 
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FIG. 12: The exit velocity versus initial phase difference 
graphs in the ODE model (|4.42p at various values of (q: (1) 
0.00707; (2) 0.00548; (3) 0.00495; (4) 0.00350; (5) 0; (6) 
-0.00350; (7) -0.00424. Here Co = -W,tpo = -0.01167, 
and e = 0.0036. 



down to zero (the integrable case), then down further 
to negative values. In this process, we closely monitor 
how the fractal structure of Fig. [S£a) changes as e de- 
creases. The result is shown in Fig. [TT] Here, the —Coo 
verse tpo graphs are plotted at six decreasing e values: 
e = 0.13665,0.036,0.0036,0,-0.0036 and -0.036. We 
see that as e decreases from 0.13665 but above zero, 
primary hill sequences and the fractal regions between 
them persist and are clearly visible in Fig. [TTT l. 2, 3). 
Indeed, we have zoomed into the sensitive regions be- 
tween primary hills in each of Figs. [TTTl. 2, 3), and 
obtained higher order structures which look very similar 
to those shown in Fig. [51 As e — > + , our key observation 
is that, the peaks of individual primary hills as well as 
the nearby fractal regions collapse to sequences of points 
on the smooth —Coo curve of the integrable system (see 
Figs. [TTT 3. 4)). Closer examination tells us that, these 
sequences of points in Fig. [TlT4) are nothing but the two 
sequences of singularity points of the integrable system 
which we plotted in Fig. [5] In other words, hill sequences 
and fractal structures in the non-integrable system bifur- 
cate from the singularity points of the integrable system! 
However, this bifurcation is one-sided: as e decreases be- 
low zero, no fractal regions appear, see Fig. ITTT 5). A fi- 
nite number of primary hills, reminiscent of primary hill 
sequences for positive e values, do exist. But the whole 



graph is smooth, and it has no fractal structures inside 
(even the spike-looking parts of the graph in Fig. [TlT5) 
turn out to be smooth upon closer examination). Fur- 
thermore, as e decreases further below zero, the number 
of primary hills keeps decreasing, and the graph becomes 
more smooth, see Fig. ITTT 6). Thus, fractal structures 
are a signature of the non-integrable system (|4.42p only 
for positive values of e, not negative values of e. 

To further substantiate our claim on fractal structures 
of the non-integrable system bifurcating from singular- 
ity points of the integrable system, we tune initial con- 
ditions so that singularity points in the integrable sys- 
tem gradually disappear in the ipo interval, and check 
if fractal structures in the non-integrable system disap- 
pear as well (for small e). Specifically, we fix the Co and 
ipo values in Eq. (|4.44[) and tune the Co value, as we 
did in Fig. [TO] The e value in Eq. (|4.42[) is taken as 
e = 0.0036, which is very small. Thus, the non-integrable 
system is weakly perturbed from the integrable one. For 
the above initial conditions, singularity points of the in- 
tegrable system have been displayed in Fig. [TO] in the 
(4>o, Co) plane. We gradually decrease the Co value. For 
each Co j we numerically compute the exit velocity ver- 
sus tpo graph in the perturbed (non-integrable) system 
(|4.42p , and compare how this graph relates to singularity 
points of the integrable system in Fig. [TO] To illustrate, 
we pick seven representative Co values, which are 0.00707, 
0.00548, 0.00495, 0.00350, 0, -0.00350 and -0.00424 in 
decreasing order. These seven Co values are marked by 
horizontal dashed lines in Fig. [TO] As we can see from 
that figure, at these seven Co values, the numbers of sin- 
gularity points in the ipo interval are 1, 3, 5, oo, oo, oo, 
and respectively. For each of these seven Co values, 
the corresponding exit velocity versus i/'o graph in the 
perturbed system (|4.42|) is shown in Fig. [TO] We no- 
tice from this figure that the numbers of primary hills 
and fractal regions near these hills at these Co values are 
equal to 1, 3, 5, oo, oo, oo, and respectively — exactly 
like singularity points in the integrable system! In par- 
ticular, when singularity points in the integrable system 
disappear, so do primary hills and fractal structures in 
the weakly perturbed non-integrable system. Further- 
more, the locations of primary hills and fractal regions 
closely follow those of singularity points of the integrable 
system. Thus, the connections between them are unmis- 
takable. Fig. [TO] together with Fig. [TTJ establishes be- 
yond doubt that primary hills and fractal structures in 
the non-integrable system (|4.42|) indeed bifurcate from 
singularity points of the integrable system. 

The bifurcation of fractal structures from singularity 
points of the integrable system indicates that near such 
points, the solutions of the perturbed system (|4.42[) are 
very sensitive to initial conditions. To shed light on why 
this sensitivity occurs, we present some numerical results 
below. First, we look at the integrable system (with 
e = 0). Taking the initial conditions as (|4.44|) , evolu- 
tions of if) versus r at the singularity point -00 = 0.98325 
(marked in Fig. [5] bottom panel) and its left and right 
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near neighbors ipo — 0.92 and 1.05 are plotted in Fig. 
[T"37 a) using the solution formula (|5.3[) . An interesting 
feature about these evolutions is that for initial ■f/'o val- 
ues at the two sides of the singularity point, the phase 
functions ip(r) have drastically different trajectories as 
they go through the time r 700 where the two soli- 
tary waves interact strongly (this time is the singular- 
ity time of the singular solution at -0o = 0.98325). For 
ipo below the singularity point, the phase sharply (but 
continuously) decreases by 2ir, while for tfj above the 
singularity point, the phase sharply (but continuously) 
increases by 2tt. Recall that ip oc A/3 and it determines 
the relative energy (amplitude) distributions among the 
two solitary waves [see Eqs. (|4UT1) . P~10")) and (|OTj) ]. 
we know that on the two sides of the singularity point, 
the energies have opposite distributions among the two 
solitary waves during their strong interactions. However, 
after the interaction is completed, the asymptotic slopes 
of the three 4>(t) trajectories in Fig. [T3T a) are almost the 
same, signaling that the interaction outcome is actually 
insensitive to the ipo values (the roughly 2tt difference be- 
tween these phase trajectories does not affect the physical 
solutions). This is why outcomes of weak interactions in 
the integrable system (|4.42p do not exhibit sensitive de- 
pendence on initial conditions (see Fig. [5] bottom panel). 
However, when system (|4.42[) is positively perturbed, the 
results are completely different. To demonstrate, we take 
e = 0.0036 now, while the initial conditions (|4.44p remain 
the same. In this slightly perturbed system, the solution 
develops finite-time singularity at the singularity point 
tpo = 0.9695, which is the counterpart of the singular- 
ity point mentioned above in the integrable system. The 
phase function at this singularity point is plotted in Fig. 
fT3lb) (solid line). (It is noted that in this perturbed case, 
we do not have exact solution formulas, hence this solu- 
tion was obtained by numerically integrating Eq. (|4.42|1 . 
Due to the finite-time singularity in the solution, our 
numerical integration can not go beyond the singularity 
time r « 700. The solution beyond the singularity time, 
shown in Fig. I13f b) as dotted lines, was inferred from our 
numerics at nearby ipo values.) On the two sides of the 
singularity point, we select two nearby values ipo — 0-92 
and 1.01. The phase trajectories at these O values are 
also plotted in Fig. [TSTb). We see that as these trajecto- 
ries go through the time r « 700, one sharply decreases 
by 2tt, while the other sharply increases by 2tt, similar 
to what happens in the integrable case (see Fig. [TUfa)). 
However, after these sharp decreases/increases, the tra- 
jectories turn around and start to move in the opposite 
direction. Eventually, these trajectories approach dras- 
tically different asymptotic slopes (one positive and the 
other one negative in fact), indicating that the interac- 
tion outcomes are very different for these slight changes 
in the V'o values. This is the phenomenon of sensitive 
dependence on initial conditions which occurs in the per- 
turbed system (|4.42[) (with e > 0), but not the integrable 
system (with e = 0). 

It is also enlightening to look at this sensitive depen- 



dence on initial conditions from the viewpoint of PDE 
evolutions. To illustrate, we take the cubic-quintic non- 
linearity (j2~7|) in Eq. (f2~Tj) , and take a = 1, /3 = 1. 
Then for the e values and initial conditions used in the 
ODE simulations of Fig. [131 and in view of the variable 
rescalings (14.401) and (14.41j) . the corresponding PDE pa- 
rameters for Fig. flUT a) (the integrable case) are 7 = 0, 
A(3 = -0.066016, AV = 0, Ax = 10, and the cor- 
responding PDE parameters for Fig. [TSTb') (the per- 
turbed case) are 7 = 0.0010, A/3 = -0.066, AV" = 0, 
Axo — 10. For these PDE parameters, the PDE evolu- 
tion results (in the form of contour plots) at three Acf>o 
values corresponding to those in the ODE simulations 
of Fig. [T3] are displayed in Fig. [T3J In the integrable 
(NLS) case (top row of Fig. [T4")) , we take the three 
A(/)q values exactly the same as those in Fig. fTBTa). i.e. 
A(j) = 0.92, 0.98325, 1.05. In this case, at the lower A^ 
value, the left solitary wave retains its higher energy at 
the collision time; at the singularity point of A(f>o, the 
two waves completely coalesce at the collision time, sig- 
naling the singularity formation in the ODE system; at 
the higher A0o value, the right solitary wave gets higher 
energy at the collision time. However, after interaction, 
the two waves always separate, and the right wave always 
gets higher energy, in all three cases. Recall that before 
interaction, the left wave has higher energy, thus we can 
call these interaction outcomes transmission. In these in- 
teractions, even though the intermediate process (espe- 
cially the collision segment) rather strongly depends on 
the initial phase difference A0Oi the interaction outcome 
is insensitive to it. These PDE evolution results com- 
pletely resemble the ODE simulations in Fig. fT3Ta). In 
the perturbed (non-integrable) case, the PDE simulation 
results are quite different from the integrable ones (as in 
the ODE simulations). In the perturbed case, we take 
the three A</> values to be 0.92,0.972 and 1.01. Notice 
that the first and third of these A0o values are exactly 
the same as those in the ODE simulations of Fig. [TSTb). 
while the middle A0o value of 0.972 is slightly different 
from the corresponding ODE value of 0.9695. This slight 
difference in the middle A(f>o value is necessary in order 
for the corresponding ODE and PDE simulations to ex- 
hibit the same behaviors, and this difference is due to the 
modeling error of the PDE evolutions by the ODE system 
(|P2"|) . At A0 O = 0.92 (see Fig. H3T1)). the interaction 
outcome is similar to the integrable ones (top row of Fig. 
[T3| in that it is also transmission. But at A</>o = 0.972 
(Fig. H3I2)), the two waves strongly coalesce, then form 
an oscillating bound state. At A<j)Q = 1.01, on the other 
hand, the two exiting waves have opposite energy distri- 
butions from Fig. I13f 1); this interaction outcome can be 
called reflection. Thus, in the perturbed case, the inter- 
action outcome is sensitive to initial conditions, which 
distinctively contrasts the integrable case. Again, these 
PDE evolution results for the perturbed case completely 
resemble the ODE simulations in Fig. [TSTb). 

The above ODE and PDE simulations corroborate the 
fact that the source of this sensitive dependence on initial 
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FIG. 14: Evolutions of \U(x,t)\ in the PDE (JUTJ with cubic- 
quintic nonlinearity (|2.7[) . corresponding to the ODE simu- 
lations of Fig. 1131 Top row: the integrable (NLS) equation 
(7 = 0); (a) A(j} = 0.92, (b) A<j5 = 0.98325, (c) A</> = 1.05. 
Bottom row: the perturbed case (with 7 = 0.0010); (1) 
A0O = 0.92; (2) A0O = 0.972; (3) A^o = 1.01. The val- 
ues of other initial parameters are given in the text. 



FIG. 13: Evolutions of ip versus r at a singularity point ipo 
and its left and right neighbor points in the dynamical system 
(|4.42j) with initial conditions (|4.44|) . (a) e = (the integrable 
case); here ipo = 0.98325 is the singularity point which is 
marked in Fig. [9] (bottom panel) ; the left and right neighbor 
points are taken as ipo — 0.92 and 1.05; (b) e = 0.0036 (the 
positively perturbed case); here ipo = 0.9695 is the singularity 
point; its left and right neighbor points are taken as ipo = 0.92 
and 1.01. 



conditions in the perturbed system lies in the finite-time 
singularities of solutions in the integrable dynamical sys- 
tem (14.421) . From the PDE point of view, the origin of 
this sensitive dependence can be traced to the coalescing 
of the two solitary waves in the integrable PDE system. 
At the moment, our understanding on this sensitive de- 
pendence and fractal structures in the perturbed system 
is still very limited. For instance, we can not yet explain 
any quantitative details inside these fractal structures, 
nor can we explain why this sensitive dependence occurs 
only for one-sided perturbations of the integrable sys- 
tem (with e > 0). These are non-trivial questions which 
merit further analysis, but they are beyond the scope of 
the present article. 

The fact of primary hills and fractal structures in the 
non-integrable system (|4.42[) bifurcating from singularity 
points of the integrable system has far reaching conse- 
quences. One important consequence is that, the main 
features of primary hill sequences shown in Figs. [Ha) 
and EJa) for PDEs and ODEs can now be analytically 



explained. For instance, the life-time formula (|3.6[) for 
primary hill sequences in the weakly perturbed system 
(|4.42[) is nothing but the analogous collision-time (sin- 
gularity time) formula (|5.16p for sequences of singularity 
points in the integrable system. To make a quantitative 
comparison between these formulas, we take the initial 
conditions (|4.44|) which was used in the PDE and ODE 
simulations of Figs. [5)Ja) and[7(a). When the time rescal- 
ing (|4.4ip is recovered, the collision-time formula (|5.16[) 
of the integrable system becomes 

0.0839i„ = 2mr + n, (6.1) 

which compares very favorably with the life-time formu- 
lae ([33]) . (|37f|) and ([446)) in direct PDE and ODE simula- 
tions. The small differences in the uj and S values between 
the analytical formula (|5.16[) and the PDE/ODE ones 
(|3.6p are caused by the not-so-small value of e = 0.13665. 
As e — + 0, these quantitative differences will vanish. Re- 
garding the locations of individual hills in the primary- 
hill sequence, they are described by the formula (|5.17|) for 
singularity-point locations of the integrable system when 
e <C 1. Note that the form of this formula is different 
from all previous ones on window sequences in solitary- 
wave collisions [15J, [lfj, |2CJ, |25j, [29( . 

For e > 0, each primary hill is paired with a sensitive 
(fractal) region at its foot (see Figs. [TT1 and IT2|) . Similar 
to primary hill sequences, the locations of these fractal 
regions are described by the same formula (|5.17[) in the 
limit e — > + . 
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The fact of primary hills and fractal structures bifur- 
cating from singularity points of the integrable system 
also explains major features of interaction results in Fig. 
[Sfa) for the exponential nonlinearity (|3.8[) . We have no- 
ticed that, unlike Fig. [T] this graph has only one infi- 
nite sequence of primary hills accumulating toward the 
left (the right sequence of Fig. [T] is absent) . This phe- 
nomenon is due to the fact that for the choices of initial 
conditions for Fig. [5ja), there is only one infinite se- 
quence of singularity points in the integrable system. To 
see it, we first calculate the f,g and e values for Fig. 
[Sfa), which are found to be 

/ = 228.8211, g = 231.91770, £ = 0.01353. (6.2) 

Thus in the scaled dynamical equation (|4.42|) , the initial 
conditions corresponding to those for Fig. E{a) are 

Co = -12.13260, Co = 0, ip„ = -0.00297. (6.3) 

Notice that ipo < 2e^°, thus according to the classifi- 
cations of singularity points in the end of the previous 
section (case 2), the integrable equation (|4.42[) with the 
above initial conditions has only one sequence of singu- 
larity points in the ipo interval, accumulating to the left 
toward ipo = + . This is in perfect agreement with the 
primary-hill sequence of Fig. Ufa) from direct PDE sim- 
ulations. 

In many of the interaction results presented in this 
paper, the exit velocity versus ipo graphs have infinite se- 
quences of primary hills (see Figs. Q] and [5] for instance); 
when zooming into the sensitive regions between primary 
hills, one gets infinite sequences of secondary hills. It 
is important to understand that these two infinite se- 
quences are pure coincidence, and are totally un-related. 
Each primary hill corresponds to a particular singularity 
point of the integrable system, thus the number of pri- 
mary hills is equal to the number of singularity points in 
the integrable system. This number can be either infinite 
or finite, depending on the choices of initial conditions. 
For instance, Figs. [T2T l. 2, 3) have 1, 3 and 5 primary 
hills, corresponding to the same numbers of singularity 
points on the top three dashed lines of Fig. [TUJ On the 
other hand, at the foot of each primary hill, there is al- 
ways an infinite sequence of secondary hills (when e > 0) . 
In other words, secondary hills always exist as an infi- 
nite, not finite, sequence. For example, if one zooms into 
each of the three sensitive regions at the foot of the three 
primary hills in Fig. [T2T 2). one always gets an infinite se- 
quence of secondary hills. Thus secondary-hill structures 
are un-related to primary-hill structures. If we zoom into 
the sensitive regions between secondary hills, we always 
get infinite sequences of tertiary hills which are very sim- 
ilar to the sequences of secondary hills both qualitatively 
and quantitatively (see Figs. E^b, c)). This process can 
continue indefinitely. Thus, our conclusion is that sensi- 
tive regions between primary hills are fractal structures 
(in the sense that portions of these structures, when am- 
plified, are the same as the strcutures themselves). But 
the whole graph with primary hills is not a fractal. 



VII. APPLICATIONS TO THE GENERALIZED 
NLS EQUATIONS WITH VARIOUS 
NONLINEARITIES 

In previous sections, we have shown that for the cubic- 
quintic and exponential nonlinearities at selected param- 
eters (a = 1,7 = 0.04, f3o = 1 for the former, and 
0o = 2.3 for the latter), weak interactions of solitary 
waves exhibit hill sequences and fractal structures for a 
wide range of initial conditions, and the reduced ODE 
model (|4.42p accurately captures these interaction dy- 
namics both qualitatively and quantitatively. In this sec- 
tion, we consider a larger question: for a given form of 
nonlinearity in the PDE (|2.1[) . can it exhibit fractal struc- 
tures? For example, with the cubic-quintic nonlinearity 
(|2.7p , for what parameters a and 7 can one possibly find 
fractal structures? This question can be answered by 
applying our previous results on the ODE model (14.421) . 
For demonstration purpose, we will do so for three forms 
of nonlinearity: cubic-quintic, exponential, and saturable 
nonlinearities. 

From the analysis of the ODE system (|4.42[) in the 
previous section, we have found that fractal structures in 
weak interactions can only occur for e > 0, not for e < 0. 
Thus, once we have obtained the functional dependence 
of £ on system parameters, it will quickly become clear 
when fractal structures can arise. The analytical expres- 
sion for £ is given in Eq. (|4.41 1) . Notice that due to the 
Vakhitov-Kolokolov stability criterion [f| [33[ , the solitary 
wave is linearly stable only when Pp > 0, i.e. £ > — 1. 
Below, we will use the £ formula in (|4.41[) to calculate e 
for general system parameters in the three nonlinearities 
mentioned above. 

First we consider the cubic-quintic nonlinearity (|2.7|) . 
When a < 0, we found that Pp is always negative, i.e. 
the solitary wave is always linearly unstable. Thus we 
only consider the a > case below. In this case, it is 
easy to see from Eqs. (|2.3j) and 12. 7j) that by a rescaling 
of variables, we can make a — 0o = 1. Thus the only 
remaining parameter for this nonlinearity is 7. Using 
the analytical formula (|2.9p for P, we can obtain the 
dependence of £ on 7, which is plotted in Fig. [TBl l). 
From this graph, we see that s > when 7 > 0, and 
£ < when 7 < 0. Thus fractal structures in this cubic- 
quintic model can appear only when 7 > 0, not when 
7 < 0. If 7 = 0, this cubic-quintic model reduces to the 
integrable NLS equation, and the dynamical equations 
(|4.42j) reduce to the integrable case (with e = 0) studied 
in Sec. 5. In this integrable case, there is of course no 
fractal dependence in solitary wave interactions. 

Next, we consider the exponential nonlinearity (|3.8p . 
In this case, the solitary wave depends only on the prop- 
agation constant [3, thus £ depends only on as well. 
The analytical expression for function e((3) is not avail- 
able, but this function can be easily determined by nu- 
merical methods, and its graph is plotted in Fig fTST b). It 
is seen that this graph has two critical propagation con- 
stants, p a = 2.2457 where £ = 0, and /?(, = 14.0051 where 
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FIG. 15: Graphs of e versus system parameters for three different nonlinearities: (1) e verse 7 for the cubic-quintic nonlinearity 
(|2.7[) : (2) e verse /3 for the exponential nonlinearity (|3.8p : (3) e verse /3 for the saturable nonlinearity (|7.1[l . 



s = +00. When /3 < f3 a , e < 0, thus fractal structures 
do not exist; when (3 a < f3 < (3b, thus fractal structures 
can appear (see previous sections); when /3 > e < — 1, 
thus the solitary wave is linearly unstable. 
Next, we consider the saturable nonlinearity, 

f(\u\ 2 ) = 1 - YT\uf i7A) 

which is common in optics (for instance, in photorefrac- 
tive crystals [HJ]). Here one is added in the above formula 
to make F(Q) = (this does not affect the solitary waves 
and their interaction dynamics). In this case, e also de- 
pends only on the propagation constant /3. This depen- 
dence is computed numerically and plotted in FigfT5FcV 
We find that e is negative for all values of (3, thus fractal 
structures can not exist in weak interactions of solitary 
waves for this saturable nonlinearity. This conclusion is 
consistent with our earlier results for the cubic-quintic 
nonlinearity, as the saturable nonlinearity (|7.1|) resem- 
bles the cubic-quintic nonlinearity (|2.7[) with a > and 
7 < 0. It is noted, however, that for the saturable nonlin- 
earity, weak interactions of solitary waves can still exhibit 
some interesting structures as shown in Figs. [Tl75, 6), 
but these structures are not fractal structures. 

From the above three examples (as well as the previ- 
ous section), we see that the reduced ODE model (I4.42[) 
enables us to accurately predict when and where frac- 
tal structures and hill sequences appear in the space of 
initial parameters of solitary waves. Based on this re- 
duced model, a global and universal understanding on 
weak interactions of solitary waves has been achieved for 
the generalized NLS equations (|2.1[) with arbitrary forms 
of nonlinearity. 

VIII. CONCLUSION AND DISCUSSION 

In this paper, we have analyzed weak interactions of 
solitary waves in the generalized nonlinear Schrodinger 
equations with general forms of nonlinearity. We have 
shown that these interactions exhibit similar fractal de- 
pendence on initial conditions for different nonlinearities. 



To analytically explain these universal fractal structures, 
we derived a set of fourth-order dynamical equations for 
the solitary-wave parameters using asymptotic methods. 
A remarkable feature of these dynamical equations is that 
they contain only one parameter, which is dependent on 
the specific form of nonlinearity. When this parameter is 
zero, these dynamical equations are integrable, and the 
exact analytical solutions are derived. When this param- 
eter is non-zero, the dynamical equations exhibit fractal 
structures which match those in the original PDEs both 
qualitatively and quantitatively. We have also investi- 
gated the origin of these fractal structures, and found 
that they bifurcate from the singularity points (i.e. ini- 
tial conditions for singularity solutions) in the integrable 
system. Based on this observation, an analytical crite- 
rion for the existence and locations of fractal structures 
is obtained. Lastly, we applied these analytical results 
to the generalized nonlinear Schrodinger equations with 
various nonlinearities such as the cubic-quintic, exponen- 
tial and saturable nonlinearities, and predictions on their 
weak interactions of solitary waves are presented. 



Regarding the bifurcation of fractal structures from the 
integrable dynamical equations, even though we have es- 
tablished that this bifurcation occurs at the singularity 
points of the integrable system, more challenging ques- 
tions are to comprehensively analyze how this bifurcation 
takes place, and to quantitatively predict the detailed 
geometric structures inside these fractals. This has not 
been done yet. Recently, Goodman and Haberman an- 
alyzed the approximate ODE models for the collisions 
of solitary waves in three physical systems where win- 
dow sequences and fractal structures have been reported 
[H, US H3| ■ They found that the origin of window se- 
quences and fractal structures in these systems lies in the 
crossing of the separatrix (homoclinic orbit). Analytical 
predictions on the locations of window sequences in the 
ODE models were derived as well. It is not clear at the 
moment whether similar analysis can be performed for 
our system (|4.42l) . This question is beyond the scope of 
the present article, and will be left for future studies. 
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